#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cutil.h>
#include "MN_cuda.h"

__global__ void MN_cuda_Kernel(int n, int p, float *vf, float *vr, float *vt, float *a, float *dm, float *m, float *kf, float *kr, float *S, float epsilon, float tolerance, float t)
{	
	if(threadIdx.x<8)	vt[threadIdx.x] = threadIdx.x;
	int tx = threadIdx.x;
	
	for(int time = 0; time< 1000; time++)
	{
		vf[tx] = kf[tx];
		vr[tx] = kr[tx];
		
		for(int j = 0; j<p; j++) {
			vf[tx] = vf[tx]*pow(m[j], (int)(S[j*n+tx]<0) ? -1*S[j*n+tx] : 0);
			vr[tx] = vr[tx]*pow(m[j], (int)(S[j*n+tx]>0) ? S[j*n+tx] : 0);
		}
		__syncthreads();
		if(tx<p)
		{
			dm[tx] = 0;
			for(int j = 0; j<n; j++) {
				dm[tx] = dm[tx] + S[tx*n+j]*(vf[j]-vr[j]);
			}
			
			m[tx] = m[tx] + a[tx]*dm[tx]*t;
			m[tx] = (m[tx]>0) ? m[tx] : 0.00001;
			a[tx] = epsilon + 1/(abs(dm[tx]) + epsilon);
		}
		__syncthreads();
	}
	vt[tx] = vf[tx]-vr[tx];
	
}


////////////////////////////////////////////////////////////////////////////////
// declarations, forward

extern "C"
//void computeGold(float*, const float*, const float*, unsigned int, unsigned int, unsigned int);


//Matrix AllocateMatrix(int height, int width, int init);
unsigned int ReadFile(float* M, char* file_name, unsigned int size);
//void WriteFile(Matrix M, char* file_name);
//void FreeDeviceMatrix(Matrix* M);
//void FreeMatrix(Matrix* M);
void steadystate(float* S, float* kf, float* kr, float* m, int n, int p);

//void MatrixMulOnDevice(const Matrix M, const Matrix N, Matrix P);

/*float S[48] = {
				1,-1,-1,0,0,0,0,0,
				0,0,2,-1,-1,0,0,0,
				0,1,0,0,-1,0,1,0,
				0,1,0,0,0,0,-1,-1,
				0,0,0,0,1,-1,0,0,
				0,0,0,1,0,-1,0,0
			};
			
float kf[8] = {1,1,1,1,1,1,1,1};
float kr[8] = {0,1,1,1,1,0,1,0};
float m[6] = {1,1,1,1,1,1};
*/
float epsilon = 0.01;
//float *epsilon = &epsilon_c;
float tolerance = 0.00001;
//float *tolerance = &tolerance_c;
float t = 0.001;
//float *t = &t_c;
//int n = 8;
//int *n = &n_c;
//int p = 6;
//int *p = &p_c;

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv) {

	float*  S;
	float*  Kf;
	float*  Kr;
	float*  m;
	//Matrix  M;
	//Matrix  N;
	//Matrix  P;

	unsigned int n,p;
	//float* vt;
	
	
	srand(52);
	
	if(argc == 6 ) 
	{
		// Allocate and read in matrices from disk
		int* params = NULL; 
		unsigned int data_read = 0;
		cutReadFilei(argv[1], &params, &data_read, true);
		if(data_read != 2){
			printf("Error reading parameter file\n");
			cutFree(params);
			return 1;
		}

		n = params[0];
		p = params[1];

		S  = (float*) malloc(n*p*sizeof(float));
		Kf = (float*) malloc(n*sizeof(float));	
		Kr = (float*) malloc(n*sizeof(float));
		m  = (float*) malloc(p*sizeof(float));

		cutFree(params);

		(void)ReadFile(S, argv[2], n*p);
		(void)ReadFile(Kf, argv[3], n);
		(void)ReadFile(Kr, argv[4], n);
		(void)ReadFile(m, argv[5], p);

		
	}
	else
	{
		
	       printf("unable reading parameter file\n");
			
	}

	// M * N on the device
    	steadystate(S, Kf, Kr, m, n, p);
   /* 
	printf("GPU computation complete\n");
    	// compute the matrix multiplication on the CPU for comparison
    	Matrix reference = AllocateMatrix(P.height, P.width, 0);
    	computeGold(reference.elements, M.elements, N.elements, M.height, M.width, N.width);
        
	printf("CPU computation complete\n");
    	// in this case check if the result is equivalent to the expected soluion
    	CUTBoolean res = cutComparefe(reference.elements, P.elements, P.height*P.width, 0.001f);
    	printf("Test %s\n", (1 == res) ? "PASSED" : "FAILED");
    
    	if(argc == 5)
    	{
		WriteFile(P, argv[4]);
	}
	else if(argc == 2)
	{
	    	WriteFile(P, argv[1]);
	}   
*/
	// Free matrices
    	//FreeMatrix(&S);
    	//FreeMatrix(&Kf);
    	//FreeMatrix(&Kr);
	//FreeMatrix(&m);
	return 0;
}

////////////////////////////////////////////////////////////////////////////////
//! Run a simple test for CUDA
////////////////////////////////////////////////////////////////////////////////
void steadystate(float* S, float* kf, float* kr, float* m, int n, int p)
{
	//Matrix P;
	float *vt;

	vt = (float*) malloc(n*sizeof(float));
	
	float *S_d;
	float *kf_d;
	float *kr_d;
	// float *epsilon_d;
	// float *tolerance_d;
	// float *t_d;
	float *vf_d;
	float *vr_d;
	float *vt_d;
	float *dm_d;
	float *m_d;
	float *a_d;
	//int *n_d;
	//int *p_d;
	
	
    	cudaMalloc((void**)&S_d, n*p*sizeof(float));
	cudaMalloc((void**)&kf_d, n*sizeof(float));
	cudaMalloc((void**)&kr_d, n*sizeof(float));
	// cudaMalloc((void**)&epsilon_d, 1*sizeof(float));
	// cudaMalloc((void**)&tolerance_d, 1*sizeof(float));
	// cudaMalloc((void**)&t_d, 1*sizeof(float));
	
	cudaMalloc((void**)&vf_d, n*sizeof(float));
	cudaMalloc((void**)&vr_d, n*sizeof(float));
	cudaMalloc((void**)&vt_d, n*sizeof(float));
	cudaMalloc((void**)&dm_d, n*sizeof(float));
	cudaMalloc((void**)&m_d, n*sizeof(float));
	cudaMalloc((void**)&a_d, p*sizeof(float));
	//cudaMalloc((void**)&n_d, 1*sizeof(int));
	//cudaMalloc((void**)&p_d, 1*sizeof(int));
	
	
	cudaMemcpy(S_d, S, n*p*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(kf_d, kf, n*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(kr_d, kr, n*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(m_d, m, n*sizeof(float), cudaMemcpyHostToDevice);
	// cudaMemcpy(epsilon_d, epsilon, 1*sizeof(float), cudaMemcpyHostToDevice);
	// cudaMemcpy(tolerance_d, tolerance, 1*sizeof(float), cudaMemcpyHostToDevice);
	// cudaMemcpy(t_d, t, 1*sizeof(float), cudaMemcpyHostToDevice);
	//cudaMemcpy(n_d, n, 1*sizeof(int), cudaMemcpyHostToDevice);
	//cudaMemcpy(p_d, p, 1*sizeof(int), cudaMemcpyHostToDevice);
	
	cudaMemset(a_d,1,p*sizeof(float));
	
	
	
	// Setup the execution configuration
	dim3 dim_grid, dim_block;
	cudaError_t cuda_ret;
	
    	// Launch the device computation threads!
		/* Kernel Call */
		cudaDeviceProp dev_prop; 
		cudaGetDeviceProperties( &dev_prop, 0);
		printf("Max Threads per block = %d\n", dev_prop.maxThreadsPerBlock);
		printf("Shared Memory per block = %d\n", dev_prop.sharedMemPerBlock);
		printf("Max threads per multiprocessor = %d\n", dev_prop.maxThreadsPerMultiProcessor);
		printf("MultiProcessor count = %d\n", dev_prop.multiProcessorCount);
		
		//int TILE_WIDTH = 16;
		//int MATRIX_WIDTH = P.width;
		//int MATRIX_HEIGHT = P.height;
		
		//printf("Tile Width = %d\n", TILE_WIDTH);
		//printf("Matrix Width = %d\n", MATRIX_WIDTH);
		//printf("Matrix Height = %d\n", MATRIX_HEIGHT);
		
		dim_block.x = n;
		dim_block.y = 1; 
		dim_block.z = 1;
    	dim_grid.x = 1;
    	//if(MATRIX_WIDTH % dim_block.x != 0) dim_grid.x++;
    	dim_grid.y = 1;
    	//if(MATRIX_HEIGHT % dim_block.y != 0) dim_grid.y++;
    	dim_grid.z = 1;

	MN_cuda_Kernel<<<dim_grid,dim_block>>>(n, p, vf_d, vr_d, vt_d, a_d, dm_d, m_d, kf_d, kr_d, S_d, epsilon, tolerance, t);

	
	/* Synchronize threads */
	cuda_ret = cudaDeviceSynchronize();
    	if(cuda_ret != cudaSuccess) printf("Unable to launch/execute kernel");
		
		
    	// Read P from the device
    	//CopyFromDeviceMatrix(P, Pd); 
		cudaMemcpy(vt, vt_d, n*sizeof(float), cudaMemcpyDeviceToHost);

    	// Free device matrices
    	// FreeDeviceMatrix(&Md);
    	// FreeDeviceMatrix(&Nd);
    	// FreeDeviceMatrix(&Pd);
		for(int i = 0; i<n; i++) {
			printf("vt[%d] = %f\n",i,vt[i]);
		}
}


unsigned int ReadFile(float* M, char* file_name, unsigned int size)
{
	cutReadFilef(file_name, &M, &size, true);
	return size;
}

